The Effects of Inter-particle Attractions on Colloidal Sedimentation 



o 
o 

c 



o 



O 

o 



o 

m 

O 

o> 

O 



x 



A. Moncho-Jorda, 1 A. A. Louis, 2 and J. T. Padding 3 

1 Departamento de Fisica Aplicada, Facultad de Ciencias, 
Universidad de Granada, Campus Fuentenueva S/N, 18071 Granada, Spain. 
2 Rudolf Peierls Centre for Theoretical Physics, 1 Keble Road, Oxford 0X1 3NP, United Kingdom 
Computational Biophysics, University of Twente, 
P.O. Box 217, 7500 AE, Enschede, The Netherlands 

We use a mesoscopic simulation technique to study the effect of short-ranged inter-particle at- 
traction on the steady-state sedimentation of colloidal suspensions. Attractions increase the average 
sedimentation velocity v s compared to the pure hard-sphere case, and for strong enough attractions, 
a non-monotonic dependence on the packing fraction <j> with a maximum velocity at intermediate 
<j> is observed. Attractions also strongly enhance hydrodynamic velocity fluctuations, which show 
a pronounced maximum size as a function of <j>. These results are linked to a complex interplay 
between hydrodynamics and the formation and break-up of transient many-particle clusters. 
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The steady state sedimentation of a set of spherical 
particles through a viscous solvent at low Reynolds num- 
ber is a classic problem in non-equilibrium physics that, 
despite its apparent simplicity, it is still far from com- 
pletely understood (HQ. A major difficulty for theories 
arises from the long-ranged (1/r) nature of the solvent 
induced hydrodynamic interactions (HI) that couple the 
motion of the particles in a complex way. Indeed, al- 
though Stokes [3J calculated the sedimentation velocity 
Vg of a single sphere in 1851, it was necessary to wait 
more than 120 years until Batchelor Q showed that the 
first correction for a hard sphere system in the dilute 
limit is given by v s /v a s = 1 — 6.550, where tj> is the parti- 
cle volume fraction. Although the sedimentation of one 
or two particles can be solved analytically [l| , the hydro- 
dynamic problem of three or more sedimenting particles 
is intrinsically chaotic [5J. Moreover, simple scaling ar- 
guments suggest that for a larger number of particles, 
the hydrodynamic velocity fluctuations around the aver- 
age Sv = v — v s should diverge with container size L as 
/(<5v) 2 ) ~ L 0]. These fluctuations have a purely hydro- 
dynamic origin and are induced by the chaotic trajectory 
of the particles, leading to velocity correlations (swirls) 
that propagate over large distances. Experiments show 
that the swirls grow with container size for smaller L (the 
unscreened regime) and then saturate for larger contain- 
ers (the screened regime) , but the nature and origins of 
the screening are still a source of controversy @, d, m, G3 • 

Colloidal sedimentation has widespread relevance in in- 
dustrial applications such as paints, coatings, ceramics, 
food, and cosmetics. Inter-particle interactions are criti- 
cal to tuning the desired dispersion properties. Neverthe- 
less, most theoretical work has focussed on purely repul- 
sive hard sphere (HS) particles, and much less is known 
about the effect of attractions on sedimentation. Batche- 
lor's extended his calculations [3] to include interactions 
beyond the HS limit 11]; for short-ranged potentials this 



can be approximated as [l[: 

v s /v° s ks 1 - [6.55 - 3.52(1 - B*)] <j> + O(0 2 ) (1) 

where the normalized second virial coefficient [l[ is de- 
fined as B% = B 2 /B2 S , with B^ s the virial coefficient 
calculated with the effective HS radius of the colloids. 
Eq. |T]) suggests that attractions should increase the sed- 
imentation velocity, while added repulsions should de- 
crease it compared to the pure HS case. Experiments on 
dilute suspensions with inter-colloid attractions (lil . [l3T ] 
or long-ranged electrostatic repulsions [14| are consistent 
with this picture. However, these theories and experi- 
ments are only relevant for very dilute suspensions. What 
happens to the average sedimentation velocity at larger 
volume fractions is not well understood, and almost noth- 
ing is known about the effect of attractions on velocity 
fluctuations at any particle density. 

In this letter we address these questions by applying 
a mesoscopic simulation technique based on stochastic 
rotation dynamics (SRD) [15j, ll6|, that can successfully 
reproduce h ydr odynamic fluctuations in steady-state sed- 
imentation [17| and has recently been shown to quanti- 
tatively describe colloidal sedimentation experiments, in- 
cluding complex non-linear effects such as Rayleigh Tay- 
lor instabilities 18]. The accuracy with which SRD was 



shown to reproduce colloidal experiments gives us con- 
fidence in the predictions of our simulations with inter- 
particle attractions. We are able to go beyond the dilute 
regime and find that, for short-ranged attractive poten- 
tials, Eq. (P) is accurate up to about <f> ~ 0.05, but breaks 
down for higher packing fractions. We measure, for the 
first time, the effect of attractions on velocity fluctua- 
tions in the unscreened regime, and find that the range 
of the hydrodynamic swirls is greatly increased by at- 
tractions, with a maximum around <fi fts 0.07 for stronger 
attractions. We link the increase in both average sedi- 
mentation velocity and in the size of the hydrodynamic 
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swirls to a complex interplay between the aggregation, 
fragmentation and sedimentation of transient clusters. 

In our simulations, the colloid-colloid interaction was 
modeled by a classic DLVO [1] potential: V cc (r) — 
VHs( r ) + V v dw( r ) + Vdh- The first term is a repul- 
sive HS like contribution of the WCA form: /3Vfrs(r) = 
10 [{a/r) 2n - (<r/r)" + 1/4] for r < 2 1 /« cr an d for 
r > 2 1 /"<7, where n — 24 and a = a cc , the colloidal 
HS radius. The second and third terms are the short- 
range van der Waals attraction and the repulsive Debye- 
Hiickel-like contribution respectively [H, Ea] • In order to 
overcome the singularity of the van der Waals contribu- 
tion at contact, we also introduced a cut-off distance (the 
so-called Stern layer) given by 8 = 0.048cr cc 0. Keep- 
ing the Debye screening length fixed at n = 8.96/cr cc , we 
varied the Hamaker constant and particle charge to ob- 
tain four different potentials with an attractive minimum 
at short inter-particle distance. The normalized second 
virial coefficients were B* 2 = -0.063, -0.507, -1.044 and 
— 1.416 respectively. Note that all these values are above 
B* 2 — — 1.5 in order to avoid fluid-fluid phase separation 
of the colloidal system fiol ]. 

Brownian fluctuations and HI are induced by includ- 
ing SRD fluid particles in the simulation. These parti- 
cles interact with each other through an efficient coarse- 
grained collision step that conserves mass, energy and 
momentum, so that the Navier Stokes equations are re- 
covered at the macroscopic level note that in the 
literature this method is also called multiple particle 
collision dynamics and has been widely applied to soft- 
matter simulations 2l|- The colloids are coupled to the 
SRD fluid through an interaction of the WCA form with 
f3V cf (r) = 10 [(cr/r) 2 ™ - (cr/r)" + 1/4] for r < 2 1 /™cr and 
for r > 2 1 /™er, where n — 6 and a — a c f — 0.465<r cc , 
and the ensuing equations of motion are updated with a 
standard Molecular Dynamics algorithm for the colloid- 
colloid and the colloid-fluid interactions, and with a 
coarse-grained SRD collision step for the fluid-fluid in- 
teractions . The colloid-fluid diameter a c f is lightly 
smaller than 0.5<r cc to avoid spurious depletion forces be- 
tween the colloids 1(| 17 1. This mesoscopic simulation 



technique has been shown to reproduce the correct low 
Reynolds (Re) hydrodynamic flow behaviour with an ef- 
fective hydrodynamic radius of a « 0.8er c /, as well as the 
correct thermal Brownian fluctuations and diffusion for 



colloidal suspensions. We refer the reader to Refs. [la ill I 
for further technical details and a justification of our SRD 
parameter choice. 

The simulations were performed by placing N c = 8 — 
819 colloids in a box of sizes L x = L y = I6a c f and L z — 
48<7 C / with periodic boundary conditions. The number 
of solvent particles was N — 40V/ r . ee /cr^ ~4-5x 10 5 , 
where Vf ree is the free volume left by the colloids. A 
gravitational external field g is applied to the colloids in 
the z direction in order to induce sedimentation. After 
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FIG. 1: Average sedimentation velocity, v s as a function of 
the volume fraction (f> = §p« for four different inter-particle 
attractions and for hard spheres. The dashed lines are Batch- 
elor's predictions |T} for the dilute limit. The dotted line is a 
prediction known to be accurate for hard spheres [23 ]. 



an initial transient time, the system reaches steady-state 
conditions, where the average sedimentation velocity v s 
is constant, the one-body particle spatial distribution is 
homogeneous, and no drift is observed. To keep the cen- 
tre of mass fixed, the downward volume flux of colloids 
is compensated by an equivalent upward volume flux of 
fluid. The simulation box sizes are small enough that 
we are still in the unscreened regime (l7j : larger simula- 
tions and possibly different boundary conditions [^, [ltj 
are necessar y in order to observe screening. The particle 
Re number Re = v s a/v, where v is the kinematic viscos- 
ity, was kept at Re « 0.08, which is small enough for 
the system to be in the correct low Re number Stokesian 
regime. Similarly, the Pcclct number Pe = v s a/ D c , with 
D c the colloid diffusion coefficient, was kept at Pe = 2.5, 
so that thermal Brownian noise is non-negligible. 

We begin our study with the effect of the inter-particle 
attractions on the average sedimentation velocity. Fig. [T] 
shows v s for various attractions and for different pack- 
ing fractions (f> = ^npa 3 , where p is the colloid num- 
ber density. Note that for these hydrodynamic effects, 
the correct radius to use is the hydrodynamic one. In 
the dilute limit, the simulation results are well described 
by Batchelor's prediction |T|) with no free parameters. 
For B2 > —0.86 the slope at \ow-</> is negative, but 
for B2 < —0.86 the attractions become strong enough 
to overcome the backflow-induced velocity reduction and 
give rise to a positive slope. These results suggest that, 
at least for short ranged attractions, the exact potential 
details are unimportant and v s /v1 is controlled by B2. 

At larger packing fractions the Batchelor prediction 
breaks down. Interestingly, for more attractive systems 
(B2 < —0.86) we observe a clear maximum in v s v.s. <f> 
that becomes more pronounced for stronger attractions, 
as predicted by theory [23j . To the best of our knowledge, 
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FIG. 2: Correlation length £|| parallel to the sedimentation 
direction for velocity swirls, as a function of packing fraction 
4> for different inter-particle attractions and for pure HS. At- 
tractions strongly enhance the size of the velocity swirls. 



this is the first time that a non-monotonous dependence 
on v s has been observed in simulations. 

Having established that attractions significantly influ- 
ence the sedimentation velocity, we now turn to the ve- 
locity fluctuations of the particles around the average, 
Sv = v—v s . Here we focus on spatial correlations in the z- 
direction, defined as C z (r) — (Sv z (Q)6v z (r)} / (5v z (0) 2 ) , 
where (...) is the average over time and over all colloids. 
The distance vector r may be taken parallel to the sed- 
imentation, C z (z), or perpendicular to it, C' z (x). The 
correlation function C' z (z) exhibits an exponential decay 
C z (z) ~ exp(— z/£||), where £| is the correlation length 
in the direction parallel to the sedimentation. It pro- 
vides a measure of the size of the hydrodynamic swirls. 
C z {x) decays as well, but also shows an anti-correlation 
region where the swirl moves in the opposite direction. 
The qualitative shape and decay of the swirls is similar to 
that observed in experimental observations and sim- 
ulations [13] for HS particles. However, as can be seen 
in Fig. [51 the correlation length £| is greatly enhanced 
by the attractions. Whereas for weaker attractions £|| 
decreases with <f>, similar to what is observed for HS par- 
ticles 0, Ell j for more strongly attractive systems the cor- 
relation length first increases with <p and then decreases, 
giving rise to a maximum around <f> w 0.07 for < — 1. 
A similar non-monotonic behavior with <fi is observed for 
C z (x) (not shown), i.e. the depth of the anticorrelation 
region grows for increasing particle attractions and, for 
the particular case of i?2 < — 1, it shows a maximum at 
intermediate colloidal density. 

Although the attractions in our system are not strong 
enough to form permanent clusters, they do enhance the 
probability for particles to cluster together in a transient 
fashion. This effect can be clearly seen Fig. [3] where we 
plot the probability distribution of transient clusters P(i) 
and their average life-time r(z) as a function of the cluster 




FIG. 3: (a) The probability of finding a transient cluster of 
size i and (b) the average cluster life-time (normalized by 
the Brownian time tb = a 2 /D co i) as a function of i for hard 
spheres and for four different inter-particle attractions. The 
particle volume fraction is cj> = 0.0233. 



size, i. In order to distinguish whether a pair of parti- 
cles belong to the same cluster or not, a cut-off distance 
r cut = l-06er cc was used, which is a reasonable estimate 
of the range of the attractive potential well. We checked 
that our results do not qualitatively depend on the exact 
cutoff distance r cut . As the attraction strength increases, 
both the probability of finding clusters (P(i > 1)) and 
the average cluster life-time t(i) increases. The grav- 
itational force on a cluster increases linearly with the 
number of particles inside the cluster, but the friction 
increases roughly linearly with the radius of gyration 
of the cluster. Since the latter typically increases less 
quickly than the former, larger clusters sediment faster 
than smaller ones, an effect we observe by tracking the 
clusters in time. Furthermore, we find that stronger in- 
teractions lead to slightly more compact clusters, with a 
smaller radius of gyration and so even larger sedimenta- 
tion velocities. This cluster picture compliments the ob- 
servation [l[ that a ttractions decrease the contribution of 
the hydrodynamic backflow term used to derive Eq. (fT|). 
Together, these effects help explain why at a fixed <j>, 
increasing the strength of the attractions increases the 
average sedimentation velocity v s . Furthermore, faster 
clusters with a longer life-time are able to move along 
larger distances with a roughly constant velocity. This 
leads to the propagation of the correlations along larger 
distances and so, to an increase of £| which tracks the 
increase in the size of the swirls. 

In Fig. [5] we show how the cluster distributions and 
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FIG. 4: The plots show (a) the probability of finding a tran- 
sient cluster of size i and (b) the average cluster life-time 
(normalized by the Brownian time tb) as a function of i for 
three different particle volume fractions, (f>. The normalized 
second virial coefficient is _BJ = —0.507. 

life-times change with 4> for a fixed attraction. Although 
we only show one value of B2 , the curves show the same 
qualitative behaviour for other values of the attraction 
we studied. For larger <fi, particles are on average closer 
to each other, and so transient clusters are more likely 
to occur, increasing v s . On the other hand, the clusters 
have a shorter life-time for larger </>, presumably because 
the collision rate with other particles increases. This ef- 
fect will decrease v s . The backflow also increases with 
4>, decreasing the average v s . At low <f> and stronger at- 
tractions, the enhanced transient cluster formation wins 
out, but at higher <j> the shorter cluster life-times as well 
as the enhanced backflow reduce v s again. In summary 
then, the competition between larger clusters and shorter 
life-times may explain the more complex dependence of 
the sedimentation velocity and velocity fluctuations on 
volume fraction that we observe for attractive particles. 

Finally, this study also raises a number of further ques- 
tions. Firstly, our simulations have been for a low Pe 
number, but it would be interesting to see what happens 
for larger Pe numbers. Preliminary simulations up to 
Pe= 15 show that, in contrast to the pure HS case where 
the v s /v® v.s. <f> curves all scale onto the same form for 
Pe numbers in this range [Tt| . we observe that v s /vq is 
slightly reduced at higher Pe number. A cluster analysis 
reveals that larger clusters are slightly less likely to occur 
at higher Pe numbers, most likely because the enhanced 
shear forces break them up, leading to a lower sedimenta- 
tion velocity. What happens in the non-Brownian infinite 



Pe number limit remains to be investigated. 

Secondly, it would be very interesting to study larger 
box sizes, to see whether attractions change the crossover 
from the screened to the unscreened regime. Thirdly, for 
even stronger attractions (-BJ < —1-5) permanent clus- 
ters should begin to form. These will then sediment more 
quickly than monomers. But as they grow and acceler- 
ate, at some point the shear forces should break them up 
again. Such a rich interplay between aggregation and hy- 
drodynamics should lead to new steady states with a clus- 
ter population that depends on the attraction strength. 
New simulations are planned to address these questions. 
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